setwd("/Users/benjaminleo/Downloads/PGMs_and_CW_onset_Empirics")


# PGMs and CCs Matching Analysis  -----------------------------------------
  #This file is for rerunning my analyses using matching techniques 

#Full data created in separate do file 
load("/Users/benjaminleo/Downloads/PGMs_and_CW_onset_Empirics795_cleaned_data.Rdata")
View(full_data)

library(tidyverse)
library(optmatch)  
library(BalanceCheck)
library(MatchIt)
library(haven)
library(margins)
library(texreg)



# Creating the Matched Data  ----------------------------------------------

match_data <- full_data
#defining relevant variables 
matching_vars <- c( "cw_onset",
                    "lag_presence",
                    "lag_count",
                    "lag_govt_pgm",
                    "lag_not_govt_pgm",
                    "lag_govt_pgm_count",
                    "lag_not_govt_pgm_count",
                    "v2x_polyarchy",
                    "v2x_polyarchy",
                    "gdppc",
                    "log_milex",
                    "log_rents",
                    "pop",
                    "rugged",
                    "total_size",
                    "peace_years",
                    "post_cw",
                    "continent")
match_data <- match_data[complete.cases(match_data[, matching_vars]), ]
rm(matching_vars)

#nearest neighbor matching 

# Nearest Neighbor Matching  ----------------------------------------------

#presence model 
matched_1 <- matchit(lag_presence ~ v2x_polyarchy +
                       gdppc +
                       log_milex +
                       log_rents +
                       rugged +
                       total_size +
                       peace_years +
                       post_cw,
                     data = match_data,
                     method = "nearest",
                     distance = "mahalanobis")
summary(matched_1)


matched_1_data <- match.data(matched_1)

match_model_1 <- glm(cw_onset ~ lag_presence,
                     data = matched_1_data,
                     weights = weights)

texreg::screenreg(match_model_1)

# Calculate predicted probabilities
pred_probs <- predict(match_model_1)

# Calculate predicted change in probability for lag_presence
pr_change1 <- margins(match_model_1, variables = "lag_presence")

# Print results
summary(pr_change1)


#govt pgms 
matched_2 <- matchit(lag_govt_pgm ~ v2x_polyarchy +
                       gdppc +
                       log_milex +
                       log_rents +
                       rugged +
                       total_size +
                       peace_years +
                       post_cw,
                     data = match_data,
                     method = "nearest",
                     distance = "mahalanobis")

summary(matched_2)

matched_2_data <- match.data(matched_2)

match_model_2 <- glm(cw_onset ~ lag_govt_pgm, 
                     data = matched_2_data,
                     weights = weights)

texreg::screenreg(match_model_2)




#Non govt PGM 
matched_3 <- matchit(lag_not_govt_pgm ~ v2x_polyarchy +
                       gdppc +
                       log_milex +
                       log_rents +
                       rugged +
                       total_size +
                       peace_years +
                       post_cw,
                     data = match_data,
                     method = "nearest",
                     distance = "mahalanobis")


summary(matched_3)

matched_3_data <- match.data(matched_3)

match_model_3 <- glm(cw_onset ~ lag_not_govt_pgm, 
                     data = matched_3_data,
                     weights = weights)

texreg::screenreg(match_model_3)

# Calculate predicted probabilities
pred_probs3 <- predict(match_model_3, type = "response")

# Calculate predicted change in probability for lag_presence
pr_change3 <- margins(match_model_3, variables = "lag_not_govt_pgm")

# Print results
summary(pr_change3)


#Creating tables from the matching analysis

match_model_list_1 <- list(
  match_model_1,
  match_model_2,
  match_model_3
)

htmlreg(match_model_list_1,
        custom.model.names = c("Model 1", "Model 2", "Model 3"),
        custom.coef.names = c("intercept",
                              "PGM Presence (t-1)",
                              "Gov't Formed PGM Presence (t-1)",
                              "Non-Gov't PGM Presence (t-1)"),
        custom.coef.map = list(
          "lag_presence" = "PGM Presence (t-1)",
          "lag_govt_pgm" = "Gov't Formed PGM Presence (t-1)",
          "lag_not_govt_pgm" = "Non-Gov't PGM Presence (t-1)"),
        custom.note = "Note: The Dependent Variable for Each Model is Civil Conflict Onset, 
        Observations are weighted by propensity scores",
        intercept = FALSE, 
        file = "pgm_cw_matching_nn.html",
        color = "black",
        caption = "",
        digits = 3
)

# Full Matching  ----------------------------------------------------------


#presence model 
matched_1b <- matchit(lag_presence ~ v2x_polyarchy +
                        gdppc +
                        log_milex +
                        log_rents +
                        rugged +
                        total_size +
                        peace_years +
                        post_cw,
                      data = match_data,
                      method = "full")
summary(matched_1b)

matched_1b_data <- match.data(matched_1b)

match_model_1b <- glm(cw_onset ~ lag_presence, 
                      data = matched_1b_data,
                      weights = weights)

texreg::screenreg(match_model_1b)

# Calculate predicted probabilities
pred_probs1b <- predict(match_model_1b, type = "response")

# Calculate predicted change in probability for lag_presence
pr_change1b <- margins(match_model_1b, variables = "lag_presence")

# Print results
summary(pr_change1b)


#govt pgms 
matched_2b <- matchit(lag_govt_pgm ~ v2x_polyarchy +
                        gdppc +
                        log_milex +
                        log_rents +
                        rugged +
                        total_size +
                        peace_years +
                        post_cw,
                      data = match_data,
                      method = "full")

summary(matched_2b)

matched_2b_data <- match.data(matched_2b)

match_model_2b <- glm(cw_onset ~ lag_govt_pgm, 
                      data = matched_2b_data,
                      weights = weights)

texreg::screenreg(match_model_2b)



#Non govt PGM 
matched_3b <- matchit(lag_not_govt_pgm ~ v2x_polyarchy +
                        gdppc +
                        log_milex +
                        log_rents +
                        rugged +
                        total_size +
                        peace_years +
                        post_cw,
                      data = match_data,
                      method = "full")


summary(matched_3b)

matched_3b_data <- match.data(matched_3b)

match_model_3b <- glm(cw_onset ~ lag_not_govt_pgm, 
                      data = matched_3b_data,
                      weights = weights)

texreg::screenreg(match_model_3b)

# Calculate predicted probabilities
pred_probs3b <- predict(match_model_3b, type = "response")

# Calculate predicted change in probability for lag_presence
pr_change3b <- margins(match_model_3b, variables = "lag_not_govt_pgm")

# Print results
summary(pr_change3b)


#Creating tables from the matching analysis

match_model_list_2 <- list(
  match_model_1b,
  match_model_2b,
  match_model_3b
)

htmlreg(match_model_list_2,
        custom.model.names = c("Model 1", "Model 2", "Model 3"),
        custom.coef.names = c("intercept",
                              "PGM Presence (t-1)",
                              "Gov't Formed PGM Presence (t-1)",
                              "Non-Gov't PGM Presence (t-1)"),
        custom.coef.map = list(
          "lag_presence" = "PGM Presence (t-1)",
          "lag_govt_pgm" = "Gov't Formed PGM Presence (t-1)",
          "lag_not_govt_pgm" = "Non-Gov't PGM Presence (t-1)"),
        custom.note = "Note: The Dependent Variable for Each Model is Civil Conflict Onset, 
        Observations are weighted by propensity scores",
        intercept = FALSE, 
        file = "pgm_cw_matching_full.html",
        color = "black",
        caption = "",
        digits = 3
)






